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An impurity solver based on a continuous-time quantum Monte Carlo method is developed for 
the Coqblin-Schrieffer model. The Monte Carlo simulation does not encounter a sign problem 
for antiferromagnetic interactions, and accurately reproduces the Kondo effect. Our algorithm 
can deal with an arbitrary number A'^ of local degrees of freedom, becomes more efficient for 
larger values of A'^, and is hence suitable for models with orbital degeneracy. The dynamical 
susceptibility and the impurity t-matrix are derived with the aid of the Fade approximation 
for various values of A', and good agreement is found with other methods and available exact 
results. We point out that the Korringa-Shiba relation needs correction for a finite value of the 
exchange interaction. 
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1. Introduction 

Strong correlations among localized and conduction 
electrons lead to the Kondo effect in impurity systems, 
and heavy fermions in a periodic lattice. Two contrast- 
ing approaches may be used to deal with lattice mod- 
els theoretically: one is to solve the model on a finite 
cluster by a method such as exact diagonalization, the 
other involves the solution of an effective impurity sys- 
tem within the framework of dynamical mean field the- 
ory (DMFT).^ The former approach is more suitable for 
low-dimensional systems, while the latter becomes exact 
in infinite-dimensional systems. For actual heavy fermion 
materials in three dimensions, the DMFT is the simplest 
approach which can take local correlations into account. 
But even within this approximate framework, perturba- 
tive calculations fail in the most interesting parameter 
range, where fierce competition arises between local and 
inter-site correlations. Therefore, numerical approaches 
are required to reliably solve the effective impurity prob- 
lem. 

In this paper, we present a new impurity solver based 
on the continuous-time quantum Monte Carlo method 
(CT-QMC) for fermion systems, which has been origi- 
nally proposed by Rubtsov^ et al. The CT-QMC eval- 
uates the infinite sum of multiple integrals in a per- 
turbation expansion by means of a Monte Carlo proce- 
dure. The CT-QMC does not require a Trotter decom- 
position in contrast with auxiliary field quantum Monte 
Carlo methods such as Hirsch-Fye. While the CT-QMC 
was first formulated as a perturbation expansion with 
respect to the Coulomb interaction, an alternative ex- 
pansion with respect to hybridization around the atomic 
limit has been developed.'' For the Anderson model, the 
latter approach does not encounter a sign problem'* and 
empirically, it is found that even more complicated mul- 
tiorbital and cluster models can be simulated without 
encountering negative weight configurations. Therefore 
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it is a powerful impurity solver for DMFT and its exten- 
sions, and has been applied to several models.'''^"* 

The CT-QMC solvers^' ^ are not suitable for very 
strong correlations. In the "segment" picture of ref. 3, the 
simulation of the Anderson model with strong Coulomb 
repulsion and a deep local level involves short excur- 
sions to high-energy states (short overlapping segments 
or anti-segments). Without special precautions, the ac- 
ceptance rate of the Monte Carlo moves will be very low. 
However, the empty and doubly occupied states play an 
essential role and give rise to the effective exchange in- 
teraction. Since these high-energy states are easily dealt 
with by ordinary perturbation theory, it is reasonable 
to employ an effective model which eliminates these vir- 
tual processes from the start. It is well known that the 
Coqblin-Schrieffer (CS) model corresponds to such a lo- 
calized limit of the Anderson model. ^ The purpose of this 
paper is to formulate the CT-QMC in a form which is 
applicable to the CS model. 

Our scheme, to be presented below, has the following 
features: 

(i) As long as interactions are antiferromagnetic, the 
scheme is free from the minus sign problem. This 
is consistent with the fact that the CS model with 
the antiferromagnetic exchange is derived from the 
Anderson model. 

(ii) Acceptance probabilities for the random walk are 
higher compared to the case of the Anderson model. 
This is because the formulation excludes the charge 
degree of freedom. The remaining magnetic ex- 
change processes are sampled efficiently. 

(iii) An arbitrary number N of local states can be dealt 
with in a simple manner. Larger values of TV make 
the computation faster in general, since the number 
of operators in each component decreases. Hence the 
scheme can easily handle orbital degeneracy. 

Since the algorithm does not encounter the minus sign 
problem, highly accurate dynamics can be derived. The 
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algorithm is therefore especially useful for large- sys- 
tems, where other methods such as exact diagonalization 
are not practical due to the rapidly increasing number of 
basis states. 

This paper is organized as follows. In the next section, 
the CT-QMC formalism is presented for the CS model. 
Section 3 discusses the Monte Carlo sampling procedure. 
The algorithm is also applied to the Kondo model which 
is discussed in §4. Numerical results for static quantities 
are shown in §5 for a wide range of N. Single-particle 
and two-particle spectra are evaluated with the aid of the 
Pade approximation in §6. The summary of the paper is 
given in §7. 

2. Formulation 

2. 1 Partition function 

In order to present our new algorithm, it is conve- 
nient to summarize the general formulation of the CT- 
QMC.^' The partition function for the Hamiltonian 
H = Hq + Hi is given by 



Z = Tl{ TrB 



-0Ho 



exp 



/ driJi(T) 
Jo 



(1) 



where Hi{t) = e^^^ifie 



-fJHo 



is an operator in the inter- 



action picture. Expanding the exponent we obtain an- 
other expression for the partition function: 



fe=o 



dn 



drfe 



X l^{Tre-^"'>Hi{Tl) ■ ■ ■ Hi{Tk)}. 



(2) 



In the CT-QMC, this sum of multiple integrals is sam- 
pled stochastically. 

We consider the Coqblin-Schrieffer model with N com- 
ponents® 



fea 



Ho = Hc + Hf 

Hi = ^ Jaa'-^aa'{ CctC^/) 



a 

(3) 



where Cfe is the energy of conduction electrons with re- 
spect to the chemical potential, and Xaa' is the X- 
operator on the local states \a) defined by Xna/ = 

— 1/2 

|Q;)(a'|. Ca = Xq X^fe'^fe"' with A^o being number of 
sites, is the annihilation operator for the conduction elec- 
tron at the impurity site. The order of the conduction 
operators in Hi specifics the definition of the equal-time 
Green function used in the perturbation expansion. In 
order to distinguish the impurity contribution to the par- 
tition function, we factorizc Z as Z = ZfZc, where Z^ is 
the partition function without the impurity. The impu- 
rity contribution is given by 

oo „/3 
fc=0 ^"^^-1 aiai a^aj, 

X S^{Trci{T'MT'{) ■ . ■ ci{Ti)c^{T'l)U 



X Trf{Tre-^"fX^,„.^ (n) • • • (rfe)}, (4) 

where conduction electron operators are grouped by 
component index, and averaged separately by (• - - )c = 
Z^^Tic{e~^^'' ■ ■ ■}■ A resultant sign in the permutation 
is represented by s. ka is the number of operators cj^Ca 
for each component a, and J^a ~ ^- Using Wick's 
theorem for the conduction electrons, we obtain 



a 

X TTf{Tre-P"tX^,^,^{Ti) ■ --X^.^^^in)}, (5) 

where /D[fc] denotes the sum of k and {ai} and the 
multiple integrals over {t^}. The ka by ka matrix D^a"^ 
is defined by {D\t°''')ij = gaij-' - rj) with Qair) = 

-{TrCaiT)cl)c. 

In the CT-QMC, statistical averages are evaluated by 
means of, for example, the Metropolis algorithm with the 
weight of Monte Carlo configurations given by eq. (5). 
The implementation of the random walk will be dis- 
cussed in §3. During the simulations, it is enough to store 
only the inverse matrix M^"^ = {D^°'^)~^. The matrix 



M, 



(fe„) 



is utilized to evaluate determinant ratios and is 



efficiently updated using fast-update formulae.^ 

2.2 Impurity t-matrix 

The conduction electron Green function Ga{T, r'; {tj}) 
for each configuration {r,} is represented as follows:^ 

Ga{T, t'; {n}) = ga{T - t') 

-E5«(^-^.)(^^^Wa(n-r'). (6) 

ij 

This equation can be derived by using the fast-update 
formula for Ma ° . An average over the Monte Carlo 
ensemble gives the physical Green function: G{t,t') = 
{G{t, T';{Ti})) MC- After the Fourier transform, we ob- 
tain 



Ga(ie„) = Sfa(ien) + 5a(ien)ta(ien)5a(ien), 



i„(ie„) = -r(^(M('=°)),,e-"(-^-^^^ 



(7) 



(8) 



MC 



where e„ = (2n + l)7rT is the fermion Matsubara fre- 
quency. Since numerical summations over i and j are 
time consuming, it is more convenient to measure in 
imaginary time as follows: 

ta{r) = -T (y2iMi'-%6ir, r, - n)\ . (9) 

\ ij I MC 

In numerical calculations, the (S-function is replaced by a 
rectangular function with a finite width, and r is sampled 
in T > with use of the anti-periodicity as proposed in 
ref. 3. 
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We should note that the t-matrix in frequency space 
inchides a constant term ta^ in the first Born approxi- 
mation 



+(1) ^ J (X ] 



(10) 



Although Ma"^ contains information of ta' at Tj = Tj, 
it is convenient to add it after the Fourier transformation 
rather than measure it in T-spacc. The constant is rele- 
vant for the calculation of the Kondo model as will be 
discussed in §4. If the f-matrix is sampled in frequency 
space, the constant term is automatically included in 
ia(ien)- 

2.3 Two-particle correlation function 

We consider two-particle correlation functions Xaa' {t) 
defined by Xaa' {t) — {T^X^^ (t)Xq/q/ ) , where the super- 
script H indicates the Heisenberg operator, and the tilde 
indicates deviation from the expectation value: Xaa = 

Xaa - {Xaa). In the CT-QMC, {TrX^^{T')X^,^,{T")) is 

evaluated by 



(1) 



2.4 Specific heat 

Thermodynamic quantities can be evaluated from the 
single-particle Green function. The expression for the in- 
ternal energy is obtained from the equation of motion for 
G(t, t') in the limit t' t + 0.^" A contribution -Eimp 
of the impurity to the internal energy is given in terms 
of the impurity t-matrix fc«(ien) for forward scattering as 
follows: 

-Eimp = {H) — {Hc)c 



Trf{Tre-0HfXa,ai (n) • • • Xa,a', (rfc)} 



(11) 



In actual computations, it is not necessary to evaluate 
the matrix products for the trace. Instead, it is suffi- 
cient to judge whether Xaa{T') and Xa'a'{T") are on 
segments of a and a' components respectively. In other 
words, we test whether Xaair') and Xa'a'ir") are per- 
mitted by the X-opcrators in front and behind of them. 
By averaging over r' with t = r' — r" , Xaa' {t) is ob- 
tained with high accuracy. The equal-time correlation 
can be represented in terms of the mean occupation by 

Xaa'(O) = {Xaa){Saa' ~ {Xa'a'))- 

We also consider response functions for M = 

'^aXaa 



X(t) = {TrM^{T)M). 



(12) 



By choosing proper m„ with m,a = 0, we can deal 
with magnetic, quadrupole and other moments within 
A'' states. Provided the local levels are degenerate and 
the system has SU(A^) symmetry, the susceptibility is 
given by the Curie constant Cm = -^"^X^a^^a- "^^^ 
dynamical susceptibility xiT~) is then given in terms of 
Xaa'{r) by 

Xir) 



N 



(13) 



a^a' 



The sums over a and a' improve the accuracy of sam- 
pling, although Xaa and Xaa' are identical for all a and 
for any combination of a ^ a', respectively. The equal- 
time correlation is x(0)/CAr = 1. The static susceptibil- 
ity is evaluated by integrating x(t). Although the static 
susceptibility can also be obtained by measuring the ex- 
pectation value of M in the presence of a small external 
field, the evaluation using x(t) gives results with higher 
precision. 



= E 



n \ 



(ie„) ia(ien)e" 



(14) 



where (5 is a positive infinitesimal quantity. 

The specific heat C is evaluated from the difference of 
-E'imp at different temperatures. Assume that the inter- 
nal energies at temperatures Tq and Ti (Tq < Ti) are 
obtained as Eq and Ex, respectively. The specific heat at 
\To + Ti)/2 is given by 

El — Eo 



7 



MC 



c 



Ti - To ' 

and its statistical error AC is estimated by 



AC = 



V(Ai?o)2 + [AE^y 



(15) 



(16) 



Ti-To 

where AEq and AEi are standard deviations of Eq and 
El, respectively. In this derivation, Gaussian distribu- 
tions have been assumed. We note that the specific heat 
computed in eq. (15) includes, in addition to the statis- 
tical errors, an error due to the finite differences, which 
is proportional to powers of Ti — Tq. Hence, in order 
to obtain a reasonable numerical accuracy, the statisti- 
cal errors of the internal energy need to be smaller than 
El — when the temperature difference Ti — Tq is de- 
creased. 

3. Monte Carlo Procedure 

In the CT-QMC method, we evaluate the statistical 
average of physical quantities by samplings with respect 
to the weight Wk in eq. (5). The random walk in con- 
figuration space must satisfy the ergodicity and detailed 
balance conditions. 

Updates which change the order of J are required for 
ergodicity, and updates which shift one of the operators 
increase sampling efficiency. In this section, we introduce 
two different algorithms for the random walk. The first, 
which changes the perturbation order by ±1 in each up- 
date, is most efficient. If some coupling constants are 0, 
however, the sampling may not be ergodic. For exam- 
ple, when the interaction lacks diagonal elements in the 
N = 2 model, the perturbation order must be changed by 
±2, because only even powers of J contribute to the par- 
tition sum. In the general case, if only some coupling con- 
stants are finite, one needs to manipulate several opera- 
tors (up to N) in one update. The algorithm is presented 
in the latter of this section. Either algorithm should be 
chosen depends on the model. 
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Fig. 1. A diagram representing the configuration of {r^} and {a^} 
of order J*. 



a„ 



a 



Fig. 2. Illustration of an insertion of a segment. 

3.1 Method 1: manipulation of a segment 

A certain configuration of order is represented in 
terms of {rj = (ti,-- - ,Tfc) and {aj = (ai,--- ,Q;fc). 
Tliese variables define tlie sequence of X-operators 

Xat,ak-i{Tk) ■ ■ ■ Xaia^-iin) ■ • • (ti ) . (17) 

Tlie corresponding c-operators are given by 
(-l)'+'4.(n)c.,(rfc)---4(^.+i)ca.(rO 

X •••4(r2)c„i(Ti). (18) 

We represent tfie above configuration by a diagram as 
shown in Fig. 1. Tfie most efficient update for {xi] and 
{ai] is the addition or removal of a single element. We 
consider the process of adding t and a, which are ran- 
domly chosen in the range [0, (3) and from the N compo- 
nents, respectively. If r satisfies t„-|_i > r > t„, {ri} 
and {ai} change into (n, • • • , r„, t, t„+i, • • • , t^) and 
(ai, • • • , a„, a, an+i, • • • , ctfe), respectively. Then one of 
the X-operators is altered as 



^Oti-i-iQti (^Ti+l ) ^ ^an + ia ('^n+1 )^aaTi ('^) 1 



(19) 



which corresponds to the change illustrated in Fig. 2. 
Namely, the segment a is inserted between a„ and an+i 
with shortening of the segment an ■ In the corresponding 
removal process, we erase one randomly chosen segment. 

According to the detailed balance condition, the ratio 
of the transition probabilities should be 

p{k^k + 1) ^ WW JV^ 

pik + l-^k) Wk k + 1' ^ ^ 

The factors N and /3 are due to the random choices of a 
and T, respectively, and fc -I- 1 due to that in the removal 
process. Since Wk has the dimension of J*^, Wk+if3/Wk is 
a dimensionless quantity. For k y^O, the ratio Wk+i/Wk 
is given by 



■exp [-?(£'„ - EaJ^ 



det 1)1+^ det Do 
det D„ det D„ 



(21) 



where I = r„+i — r is the length of the new segment. 
is the matrix with c]j(t„+i)cc((t) added to Da, and Da„ 
is the matrix with one of the operators shifted in time 
according to cJj,^(t„+i) cJj^(t). The ratio of determi- 
nants can be evaluated using fast-update formulae.^ If 
a = a„ in Fig. 2, the change is just an addition of a 
diagonal element Xaa{T), so that eq. (21) is reduced to 

(+) 

(22) 



Wk^ 



1 _ _ det 



Wk detD„ ■ 

Here d'^^ is a matrix in which cJ,(t)cq(t-(-0) is added to 
the original one. The equal-time Green function in 
should be ga{+0) to keep the probability positive. This 
is verified by taking the CS limit in the corresponding 
formula of the Anderson model (see Appendix A). In the 
case of fc = 0, we should average over the local states 
since all states contribute to the trace. Then Wi/Wq is 
given by 

Wi 

= -JaaPaga{+0), (23) 



Wn 



by 



where pa is defined 

exp{-pEa)/j:a' eM-PEa'). 

The ratios of the weights in eqs. (21)-(23) change 
their signs depending on the signs of the coupling con- 
stants. We have confirmed by numerical calculations that 
the probability remains positive in the case of antiferro- 
couplings, i.e., Jaa' > 0. This is consistent with the fact 
that the CS model with antiferro-couplings is derived 
from the Anderson model, where the minus sign problem 
does not appear. We also note that we could consider an- 
other process in which a is inserted between q;„_i and 
a„. It corresponds to a diagram where a and are ex- 
changed in Fig. 2. However, the ergodicity condition can 
be satisfied without this process. 

3.2 Method 2: manipulation of a set of operators 

If some coupling constants are zero, the ergodicity may 
require updates that change the diagrams by several per- 
turbation orders. At most N operators should be inserted 
and removed in one update. To this end, we define irre- 
ducible sets of operators, which consist of k operators 
with a set of k distinct components (a'^, • ■ • , a',,). Then 
the irreducible sets of AT-operators are given by 



(24) 



We consider a process of insertion and removal of the 
irreducible operator, which is illustrated in Fig. 3. 

To determine an irreducible operator, we begin by ran- 
domly choosing k from 1 to iV and r{ from the range 
[0,f3). If t{ satisfies t„+i > t{ > Tn, is fixed as 
a'f- — an. The remaining components a'j (1 < j < k) 
are determined so as to be different from each other. 
The imaginary times (1 < j < k) are randomly cho- 
sen from the restricted range (rj_;^ , r„+i). In the case of 
a removal process, we choose one, say a„, from all the 
segments. If two identical components appear between 
the selected segment and the next segment with compo- 
nent a„ , the operator set is not irreducible and hence the 
update is rejected at this stage. If this is not the case, we 
proceed to attempt a removal of the operators. 
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a„ 



H — -, — I — ^ 



(=a„) 



a'l a„ 



Fig. 3. Illustration of an addition of an irreducible operator. 



From the detailed balance condition, we obtain the 
update probability of the above process ior k ^ as 



1 " 

— IT' 



(max) 



p{k + K^k) Wk {N-ny.k + KjJ^ 

-1 1 1 inax I r) j 7 (max) 

IS given by /j — P and ^ — 



, (25) 



where is given by l^™^^") — p and <j ' = Tn+i 

Tj_i for j > 1. The factor A^!/(A^—k)! is the inverse of the 
probability that the corresponding irreducible operator is 
chosen. The ratio Wk+n/Wk is given by 



detD^t^ 

1 ^ 

' det Da' 

j 

(26) 



where Zj is a length of the segment a'y D^^^ is a matrix to 

which cjj, (rj^]^)cc(^. (rj) has been added, with t^+i = r{. 

In the case of fc = 0, the above procedure is applicable 
with Tn ~ and r„+i — f3. One then has to choose a'^, 
at random, and take the trace over the local states. The 
transition probability for k — becomes 



p{0 k) N -Nl 



(27) 



The additional factor N originates in the random choice 
of a'l^. We note that the factor l/(fc + k) in eq. (26) is 
not necessary for fc = because no choice is made in the 
removal process. The ratio W^/Wq is given by 

Wo = "^/o n l^-;-;-. exp(-Z,£;„, ) deti^W} , 

(28) 

where we define Zjq by Zfo = exp(— /3i?Q,). 

4. The Kondo Model 

When we deal with 5* = 1/2 systems without orbital 
degeneracy, the Kondo model is more favorable than the 
CS model, because of the particle-hole symmetry in the 
former. The algorithm for the CS model is applicable to 
the Kondo model after some modifications. 

The Kondo model is given by 



H = ^ ekci^Cka + JS ■ (Tc 



(29) 



where cTc is the Pauli matrix for conduction electrons. 
The Hamiltonian is rewritten in terms of the CS inter- 
action as follows: 



where v = —J/2. Therefore the algorithm for the CS 
model is applicable by including the potential scattering 
term into Hq. 

We introduce the Green function including the poten- 
tial scattering, g{z), which is a scalar in the spin indices. 
In terms of the bare Green function g{z), g{z) is ex- 
pressed as 

9 



9 = 



1 



vg 



(31) 



In the simulation of the CS model, 17(2:) is replaced by 
g{z). Within the CT-QMC, the impurity t-matrix tj{z) 
is computed with respect to g{z). To obtain the i-matrix 
t{z) of the Kondo model, eq. (29), the contribution of 
the potential scattering should be subtracted from g{z). 
The full Green function G{z) can be expressed as 

G = 5 + gtjg = 9 + gtg. (32) 

Solving eq. (32) with respect to t{z), we obtain 

1-vg (1 - vgy 

where the first term is the t-matrix due to the poten- 
tial scattering. Concerning the two-particle correlation 
function, no modification is required, because it does not 
depend on the selection of the basis set in the perturba- 
tion expansion. 

5. Imaginary-Time Data and Static Quantities 

We apply our algorithm to a model with SU(Af) sym- 
metry. In this case, both types of updates introduced in 
§3 are available. We have confirmed that the results of 
the two methods agree within error bars. Therefore we 
use the more efficient rank one updates for the rest of 
the calculations. 

We use a constant density of states for the conduction 
electrons 

--Im.g„(e + i(5) = po^(^-|e|), (34) 

where 6{e) is a step function and po ~ 1/2D. In the 
Matsubara formalism, the Green function is represented 
as g(ie„) = -2ipo arctan(i:)/e„) and A^tT^ 1]^ 5fc(ien) = 
-(e^ -I- D'^)~^. We choose the unit I? = 1 in this pa- 
per. The standard deviations in the MC ensembles are 
evaluated from 20 averages. In MC simulations, we have 
observed minus sign probabilities at temperatures lower 
than the Kondo temperature. However, they appear only 
at the rate of about one in lO'' updates, and may be due 
to rounding errors. 

In the following, wc show numerical results for physi- 
cal quantities without analytic continuation, such as the 
static susceptibility and the specific heat. These results 
are, within error bars, exact. Spectral function in real 
frequencies will be shown in the next section. 
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Fig. 4. (a) Probability distribution P{k) for terms of order j'', 
and (b) t-matrix t(ien) for N = 1. The temperature is chosen as 
T = 0.01. 




Fig. 5. Dynamical susceptibility x{t) ™ the imaginary-time do- 
main for AT = 8 and J = 0.075. 



CT-QMC 
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5.1 N = 1: potential scattering 

The Hamiltonian, eq. (3), is trivially solvable when 
iV = 1, which corresponds to potential scattering. The 
exact t-matrix for the potential scattering is given by the 
first term on the right-hand side of eq. (33) with v = J, 
which takes into account an infinite sequence of scat- 
tering events. We apply our algorithm to the potential 
scattering to compare with the exact solution, prior to 
applications to > 1. Figure 4(a) shows the probabil- 
ity of appearance of perturbation terms of order j'^ in 
the Monte Carlo simulations. Although the exact results 
include infinite scries of scattering, only terms of finite 
order of J are significant in Monte Carlo simulations. 
The center of the distribution increases as J increases. 
Results for the i-matrix in the Monte Carlo ensemble av- 
erage are shown in Fig. 4(b) together with the analytical 
result. The term t^^^ given by eq. (10) is subtracted in 
this figure. The exact results arc completely reproduced, 
which demonstrates the convergence of the perturbation 
series, and shows that sampling finite perturbation or- 
ders is sufficient. 

5.2 N = 8: large-N case 

The present algorithm can handle arbitrary A^. We first 
take N = 8 and give exemplary results of static quanti- 
ties as well as raw data in the imaginary-time domain. 
Figure 5 shows the dynamical susceptibility x(t) for sev- 
eral values of T. We have plotted with lines because the 



Fig. 6. Temperature dependence of the static susceptibility x for 
N = 8 and several values of J. Dashed lines are results computed 
in the NCA. 

intervals between data points are fine enough, and the 
associated errors are negligible. It turns out that the re- 
duction of the correlation, which starts at x(r = 0) = 1, 
becomes more rapid for lower temperatures. This im- 
plies Kondo screening at low temperatures. The inset 
shows x(r) against t//3. At T = 0.001, the correlation 
almost disappears at r = /3/2, while correlations remain 
for T = 0.005. This indicates complete screening of the 
local moment at T = 0.001. 

The static susceptibility x('^ = 0) can be evaluated 
by integrating xi''')- Figure 6 shows the temperature de- 
pendence of the static susceptibility for several values of 
J fov N = 8. For comparison, we plot results computed 
in the non-crossing approximation (NCA), which incor- 
porates terms up to the next-leading order in the 1/A'^ 
expansion by integral equations. The CT-QMC re- 
sults are in excellent agreement with the NCA at high 
temperatures. On the other hand, deviations are visible 
at lower temperatures. This is due to the inaccuracies of 
the NCA at low temperatures and low frequencies.^^ We 
have confirmed in fact that the incorrect upturn appears 
in the NCA results for smaller N, where the accuracy of 
the 1/N expansion diminishes. On the other hand, the 
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Fig. 7. The impurity t-matrix ta{T) in the imaginary-time do- 
main for AT = 8 and J = 0.075. 

CT-QMC produces proper values of the static suscepti- 
bility at all temperatures. 

We next show results for the impurity t-matrix. Fig- 
ure 7 shows the f-matrix f(r) in the imaginary-time do- 
main for the same parameters as in Fig. 5. Statistical 
errors are so small that we have plotted with lines as 
in the case of xi^)- As temperature decreases, — ^(t) in- 
creases. Correspondingly, the discontinuity at the bound- 
ary t{—0) — t{+0) increases as temperature decreases. 
Since the discontinuity is given by — tt^^ J lmt{'jj)doj, the 
increase indicates enhanced scattering due to the impu- 
rity. This corresponds to the formation of the Kondo sin- 
glet. 

The impurity t-matrix with imaginary frequency de- 
scribes thermodynamic quantities as shown in eq. (14). 
The temperature dependence of the internal energy E{T) 
is shown in Fig. 8(a). The statistical errors are invisible 
in this figure. The data give the specific heat C(T) via 
eq. (15). Figure 8(b) shows the temperature dependence 
of C(r) plotted together with results from the NCA." 
Statistical errors become larger at low temperatures due 
to the small mesh of temperatures. The CT-QMC agrees 
with the NCA within error bars, and reproduces a peak 
due to the Kondo effect. 

5.3 Comparison between different N 

So far we have presented imaginary-time data and re- 
sultant static quantities for iV = 8. Next we examine 
different values of N. Coupling constants are chosen so 
as to fix NJ to yield almost equal values of the Kondo 
temperature. We define the Kondo temperature Tk by 
means of the static susceptibility x by Tk = Cn/x 
low enough temperature. Table I shows x at T = 0.001 
and calculated values of Tk for = 2, 4, 6 and 8. For all 
N, Tk is less than 1 /30 of D, and therefore can be consid- 
ered small compared to D. Figure 9 shows probabilities 
P{ko,) of appearance of for different A at T = 0.001. 
The size of the matrix is equal to /c^. It turns out that 
the peak of P{ka) shifts to higher values of k^ for smaller 
N. This is because the power of J is divided into N com- 
ponents, resulting in the decrease of ka- Consequently, a 
larger value of N reduces the computational burden, and 
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Fig. 8. Temperature dependence of the internal energy E{T) and 
tlic specific heat C(T) for = 8 and J = 0.075. The dashed Une 
is the result computed in the NCA. 



Table 1. The static susceptibility x and the Kondo temperatures 
Tk estimated by Tk = Cjv/x at T = 0.001. Values in bracket 
are the standard deviations. 



N 


J 


X(T = 0.001)/Civ 


Tk 


2 


0.300 


30.95 (0.24) 


0.0323 (0.0002) 


4 


0.150 


31.67 (0.14) 


0.0316 (0.0001) 


6 


0.100 


31.03 (0.11) 


0.0322 (0.0001) 


8 


0.075 


30.53 (0.10) 


0.0328 (0.0001) 



makes it possible to reach temperatures much lower than 
Tk. 

For static quantities, exact solutions have been ob- 
tained based on the Bethe ansatz.^^ Temperature depen- 
dences of all physical quantities are given through a sin- 
gle energy scale. Hence wc can compare our results with 
the exact solutions, using the Tk determined above. We 
note that the characteristic energy Tq in ref. 15 relates 
to our Kondo temperature Tk by Tk = 2ttTq/N. Fig- 
ure 10(a) shows the static susceptibility as a function of 
T/Tk for the parameters listed in Table. I. Dashed lines 
are the Bethe ansatz solutions for A = 2 and 8. Compar- 
ison between our results and the exact solution reveals 
that Tk determined by Tk = Cn/x systematically de- 
viates from the exact results. This is ascribed to the fi- 
nite cutoff of the conduction band in our model, which is 
shown to enhance the static susceptibility. Consequently, 
Tk determined by x tends to be smaller than the Bethe 
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Fig. 9. Probability distributions P(A:a) at T = 0.001. Parameters 
and Tk axe listed in Table I. 
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Fig. 10. Temperature dependence of (a) the static susceptibility 
and (b) the specific heat for AT = 2, 4, 6 and 8. Parameters 
and Tk axe listed in Table I. Dashed lines axe the Bethe axisatz 
solution. 

ansatz value. The effect of the finite band width will 

be discussed in detail later. Temperature dependences 
of the specific heat are shown in Fig. 10(b). Large peaks 
around T/Tk ^ 0.3 are due to the Kondo effect, while 
small peaks at about T/Tk ^ 10 originate from the cutoff 
of the conduction band. We recognize systematic devia- 
tions between our results and the Bethe ansatz solutions, 
corresponding to those found in the static susceptibility. 
However, the overall behavior agrees well if we scale tem- 



n 

Fig. 11. t-matrix t(ien) of the Kondo model and the CS model 
for J = 0.3 and T = 0.001. 

perature independently of the static susceptibility. 

5.4 The Kondo model 

We can deal with the Kondo model in the way pre- 
sented in §4. In applying our algorithm for the CS model 
to the Kondo model, we use unperturbed Green func- 
tions without particle-hole symmetry, although the orig- 
inal model is particle-hole symmetric. Therefore it is a 
strict check of our algorithm whether the particle-hole 
symmetry is recovered after solving the CS interactions. 
Figure 11 shows the t-matrix t{ien) for J = 0.3 and 
T = 0.001. For comparison, a result for the CS model is 
plotted for the same coupling constant and temperature. 
In the imaginary-time representations. Re ta(ie„) = 
indicates symmetry with respect to particle-hole excita- 
tions. It is verified that the result for the Kondo model 
keeps the particle-hole symmetry, while the CS model 
does not due to the potential scattering. 

6. Dynamical Quantities with Analytic Contin- 
uation 

We proceed to spectral functions in real frequencies. 
To perform analytic continuations, we employ two kinds 
of conventional methods: the Fade approximation^*' and 
the maximum entropy method (MEM).^^ The Fade ap- 
proximation fits data at the Matsubara frequencies in 
the upper-half plane with use of a rational function, and 
extrapolates onto the real axis. Since the Fourier trans- 
form drives away noises in imaginary-time data to high 
frequencies, reasonable accuracy is expected at low fre- 
quencies. However the extrapolation is sensitive to sta- 
tistical and numerical errors, especially at high energies 
which are far from the imaginary axis. If the data include 
significant errors, the Fade approximation does not work 
in general. In this case, the statistical errors can be taken 
into account by the MEM. The MEM derives spectral 
functions from data in the imaginary-time representa- 
tion without Fourier transformation. In this subsection, 
we shall show results by the Fade approximation, since 
our data have been obtained with sufficient accuracy. 
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Fig. 12. The t-matrix — Imt{t.i;+i<5) computed in the Pade approx- 
imation at T = 0.001 and comparison with the NRG (A^ = 2) 
and the NCA (N = 6 and 8). Parameters are listed in Table I. 
Marks at oj = indicate exact values obtained from the Friedel 
sum rule, eq. (35). 



6.1 Impurity t-matrix 

Figure 12 shows —hnt^Lu) computed from the Pade ap- 
proximation at r = 0.001 for the parameters hsted in 
Table I. We plot N components separately to check the 
accuracy of the analytic continuations. We find an excel- 
lent agreement between all components around the Fermi 



s 




Fig. 13. The t-matrix —Imt(a)-|-i(5) of the Kondo model computed 
in the Pade approximation. Parameters are J = 0.3 and T = 
0.001. The dashed line is the NRG result at T = 0. 



level, while at higher energies the spectra show some dif- 
ferences. The extent of the difference can be regarded as 
a measure of the error in the approximation. For com- 
parison, we plot results computed using the numerical 
renormalization group (NRG)^*'^'' at T = for A/' = 2, 
and results from the NCA for N = 6 and 8. It turns out 
that for all N the overall shapes obtained by the different 
methods are in good agreement. Results for the Kondo 
model arc shown in Fig. 13. In computing the spectrum, 
we have set the real part of f(ie„) to neglecting tiny 
statistical errors. We have confirmed that the spectral 
function is symmetic with respect to w = and agrees 
with the result from NRG. 

At T = the Friedel sum rule relates the f-matrix at 
the Fermi level with the occupation number of the local 
state. For the constant density of states represented in 
eq. (34), -lmt„(0 + iS) is given by^° 

—lmta{Q + i5) = sin^(7mo,), (35) 
Trpo 

where n„ = (A„q,) is the mean occupation of the state 
a. The exact values of —lm.t{0 + i6) have been marked in 
Figs. 12 and 13. We remark that our results computed in 
the Pade approximation agree with the Friedel sum rule 
for all N. We conclude that analytic continuation by the 
Pade approximation works for the CT-QMC data. In the 
case where there is an energy gap or several excitations 
in the spectrum, however, the Pade approximation may 
not provide this level of accuracy. 

6.2 Dynamical susceptibility 

We next present spectra of two-particle response func- 
tions. Figure 14 shows Imx(a; + i5)/LU computed in 
the Pade approximation for the same parameters as in 
Fig. 12. Since a comparison between N components, 
which is made in Fig. 12, is not available in this case, 
we have confirmed reproducibility of the spectra between 
different Monte Carlo ensembles. Spectra fov N = 2 are 
almost Lorcntzian for all plotted temperatures, while for 
N > 4 a. peak appears at a finite energy at low temper- 
atures. This is consistent with the t-matrix in Fig. 12, 
where the Kondo resonance shifts to higher energy with 
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Fig. 14. Imxi^)/^ for several temperatures. Parameters are 
listed in Table I. Marks at w = indicate values deduced from the 
Korringa-Shiba relation, eq. (36), with a = 1 and x{T = 0.001). 
For iV = 6 and 8, results from the NCA are plotted by a dashed 
line. 



increasing N. Wc also plot, for reference, results com- 
puted in the NCA for A'' = 6 and 8. In Figs. 14(c) and 
(d), we can see excellent agreement between two results 
at T = 0.01. which is of the order of Tk. At lower tem- 
peratures, on the other hand, both results do not agree 
around a; = due to the unphysical peak of the NCA.^^ 
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Fig. 15. The parameter a in eq. (36) as a function of NJpo for 
N = 2 and 8. 



In the Fermi liquid regime, the Korringa-Shiba relation 
(KSR) connects the imaginary part of the dynamical sus- 
ceptibility with the static ones.^"'^^ It has been proven 
for the orbitally degenerate Anderson model in the wide 
band limit. In models with finite band width, however, 
the static susceptibility is enhanced and does not satisfy 
the KSR. This is due to the fact that the cutoff of the con- 
duction band enhances the static susceptibility over the 
universal value, while it does not affect the low-energy 
imaginary part. Appendix B demonstrates the deviation 
for the non-interacting Anderson model. To account for 
deviations from the wide band limit, we introduce a pa- 
rameter a and modify the original KSR as 

lmx{u) + id) n 



lim 



(36) 



^ NCn 

Here a = 1 corresponds to the KSR, but a is greater than 
unity in the case of finite band width. We have marked, 
in Fig. 14, values of Imx(u; -I- i5)/uj\^^o deduced from x 
at T = 0.001 and a = 1 in eq. (36). Comparing with 
spectra at T = 0.001, we observe large deviations for 
smaU N. For N = 2, x with a = 1 leads to about 50% 
deviation from the actual value of Imx(w -|- i6)/u)\u,=o. 
Since the width of the Kondo resonance is not negligible 
as compared with D especially for = 2, as shown in 
Fig. 12, the effect of the finite cutoff cannot be neglected. 

To see the validity of the KSR in the wide band limit, 
or cquivalcntly, in the weak coupling limit, we plot the 
parameter a against NJpo in Fig. 15. In this figure, error 
bars have been estimated from the statistical errors of %• 
We note that the errors actually come from inaccuracies 
of Imx(w + iS) as well due to the analytic continuation. 
For N = 8, the parameter a deviates from unity only by 
about 5%. For A" = 2, on the other hand, a varies lin- 
early against NJpQ. The deviation from unity is almost 
proportional to NJpo in the range of our Monte Carlo 
simulations. Consequently, in order to satisfy the KSR, 
NJpQ should be much smaller than unity for N = 2. 

7. Summary 

We have extended the CT-QMC method for fermions 
to the Coqblin-Schrieffer model. An arbitrary number N 
of local states can be treated by our algorithm and for 
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antiferromagnetic interactions, the scheme is free from 
the minus sign problem. Therefore it is a powerful tool 
to study heavy fermion systems within the framework of 
DMFT. 

We have computed spectral functions using the Pade 
approximation for the analytic continuation. The impu- 
rity t-matrix shows an excellent correspondence between 
the N components at low frequencies, which demon- 
strates the good quality of the imaginary-time data ob- 
tained using the new scheme. The accuracy of the mag- 
netic spectra has been examined by comparison with the 
NCA for large N. Both results are in good agreement at 
temperatures of the order of Tk- At lower temperatures, 
where the NCA cannot accurately reproduce the low- 
energy excitations, the present scheme still works prop- 
erly. 

The spectra have been tested using available Fermi 
liquid relations. We have found deviations from the 
Korringa-Shiba relation, whereas the i-matrix satisfies 
the Friedel sum rule for all N. The deviations are due 
to the finite cutoff of the conduction band, which affects 
the static susceptibility but not the low-energy spectrum. 
These results reveal a weak point of the CT-QMC: it is 
difficult to approach the universal regime with regard to 
the static susceptibility. 

The present algorithm can easily be applied to DMFT 
simulations of the Coqblin-Schrieffer lattice model. 
Namely, the effective medium is optimized with use of 
the impurity t-matrix. It is possible to address the for- 
mation of heavy quasi-particles through the conduction 
electron Green function as well as the t-matrix of the 
lattice systems. 
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Appendix A: Derivation of update probabilities 
in the CS model from the Ander- 
son model 

The CS model is derived from the Anderson model 

by taking the localized limit with strong correlations.^ 
Namely, charge fluctuations are suppressed by infinite 
Coulomb repulsion and the deep local level. We take the 



limits e / 



-oo and V 



oc with J 



-V^/ef fixed. 



where e / and V denote the energy of the local level and 
strength of hybridization, respectively. Similarly, update 
probabilities in the CS model can be derived from the 
corresponding expressions in the Anderson model. In this 
appendix, we demonstrate the way to take the limit in 
the Monte Carlo formalism. 

We consider the Anderson model of spin-less fermions 
for simplicity. In this case, the localized limit leads to a 
potential scattering, or equivalently, the CS model with 
= 1. The update probability for cutting a segment 
(addition of an anti-segment) is given by^ 



where I denotes a length of a segment which will be re- 
moved. D^^^ is the matrix obtained by adding the opera- 
tors (t)c(t -I- /) to the end of the matrix D. In the limit 
of e/ ^ — oo, the probability of such an update with I 
finite becomes due to the factor e'"^^. Hence I should 
be reduced as the inverse of |e/| in the update. For this 
purpose, we introduce a cutoff length Iq defined by 



— , (A » 1) 



(A-2) 



If Z > Zo) the update is negligible due to the factor e~^. 
Hence we can restrict the length tol < Iq- The restriction 
for I replaces the factor Zmax with Iq in eq. (A-1). Taking 
I = xlo with < a; < 1, the update probability is given 
in terms of J by 

det \ 



pik^k + 1) 



JXe 



-Ax 



/3 



(A-3) 



p{k+l^k) ""^ V detD J k + 1' 

Since the probability is independent of x in the limit of 
e/ ^ — oo, we integrate out x as follows: 

nl 

da;Ae-^^ ~ 1. (A-4) 



/ 

Jo 



This equality becomes exact in the limit A — » oo, which is 
realized in the limit e/ — oc . As a result , we obtain the 
update probability that the localized electron is removed 
for an infinitesimal time as follows: 



p{k^k + l) _ J / det£>(+) \ /? 



p{k + l^k) "\ detD Jk + l ^^'^^ 

This formula is identical with the probability that the 
operator Jaa^aa(— c^c^) is added. It is obvious from 
this derivation that is the matrix with c^t)c{t + 

0) added to the original one, and that the equal-time 
Green function should be g{+0). In a similar manner, 
all formulae of transition probabilities in the CS model 
can be derived from the corresponding processes in the 
Anderson model. 

Appendix B: The Korringa-Shiba relation in a 
model with finite band width 

The Korringa-Shiba relation connects Imx(w -|- 

i(5)/w|cj=o to X^(0) by a imivcrsal value. The equation 
has been proven in the wide band limit. Hence it may not 
be satisfied in numerical calculations for systems with fi- 
nite band width, as in the present study. In order to clar- 
ify the deviation from the universal value, we consider 
the non-interacting Anderson model with A''-fold degen- 
eracy. Assuming a constant density of states, eq. (34), 
the Matsubara Green function is given by 



ef + (2iA/7r) arctan(L)/e„) 
-1 -e/-hiAsgn(e„), (B-1) 



where A = irV'^po and 



2A 

ttD' 



(B-2) 



p{k^k + 1) 
p{k + l^ k) 



det Z)(+^ \ have assumed A/D <^ 1. The wide band limit corre- 

. p. ) TTTT' (A'l) sponds to a = 1, and a finite band produces a correction 

proportional to A/D. 
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We evaluate the dynamical susceptibility using the 
above Green function including the effect of the finite 
band width. The dynamical susceptibility is given by 

Xiiun) = -iVCjvT^G/(ie„OG/(ie„' +iJ^„), (B-3) 

n' 

where Vn = ^mrT is the boson Matsubara frequency. 
Evaluating the imaginary part at T = 0, we obtain the 
well-known relation Imx(cj + iS)/uj\Lu=o — nNCiyp'fiO)- 
It turns out that the quantity a docs not affect the low- 
energy spectrum. On the other hand, the real part is 
influenced by the cutoff of the band. The sum over the 
Matsubara frequency is replaced by an integral along the 
imaginary-frequency axis, and evaluates to 

x{0) = aNCNPfi.O). (B-4) 

Therefore for o = 1, the Korringa-Shiba relation is sat- 
isfied. On the other hand, if a > 1, the static suscepti- 
bility is enhanced, and the Korringa-Shiba relation does 
not hold. In the Anderson model with ?7 7^ or in the 
Coqblin-Schrieffer model, A is replaced by the width of 
the Kondo resonance, and therefore is of the order of the 
Kondo temperature. 
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